#include "header.h"

// 用基本的kernal
__global__ void kernal_basic(unsigned int* col_index, float* val, float* x, float* y, unsigned int* WLB_row_offset, unsigned int* combine_meta_vec, int bit_flag_bit_num, int reduce_offset_bit_num, int TLB_relative_row_index_bit_num, int WLB_num, int TLB_num, int nnz_of_TLB)
{
    // 查看当前warp号，然后是warp内的线程号
    int global_tid = blockDim.x * blockIdx.x + threadIdx.x;
    int warp_id = global_tid / 32;
    // thread的线程内号
    int tid_in_warp = global_tid % 32;
    int warp_num = blockDim.x * gridDim.x / 32;

    for (unsigned int WLB_id = warp_id; WLB_id < WLB_num; WLB_id = WLB_id + warp_num)
    {
        unsigned int WLB_first_row = WLB_row_offset[WLB_id];
        
        // 执行复杂的归约过程
        // 获取合并之后元数据，从这里计算出总TLB号，每个WLB有32个TLB。
        unsigned int global_TLB_id = WLB_id * 32 + tid_in_warp;

        // 获得每一个TLB的元数据
        unsigned int combine_meta = combine_meta_vec[global_TLB_id];

        // 获取每个TLB第一个加和起始位置的相对行索引以及归约偏移，这两个数据都是先右移再左移的方式带来，也可以直接用与运算来处理。
        unsigned int TLB_first_reduce_row = combine_meta << (32 - TLB_relative_row_index_bit_num) >> ((32 - TLB_relative_row_index_bit_num));
        
        // TLB向全局内存写中间结果的行偏移量
        unsigned int y_offset = TLB_first_reduce_row;

        // 获取当前线程的归约偏移量
        unsigned int reduce_offset = combine_meta << (32 - TLB_relative_row_index_bit_num - reduce_offset_bit_num) >> (32 - reduce_offset_bit_num);

        // 用一个变量记录之前是否出现过sum起始点
        bool through_sum_begin_bit = false;

        // 用一个变量存储行脑袋和行身体的中间结果
        float row_head_tmp_result = 0;
        // 除了头之外其他部分的中间结果
        float row_other_tmp_result = 0;
        // 加和的中间结果
        float sum_tmp = 0;

        // 当前非零元的全局索引
        // assert(global_TLB_id < 160);
        unsigned int global_nz_index = WLB_id * nnz_of_TLB * 32 + tid_in_warp;
        
        // 遍历线程所有的非零元
        for (unsigned int TLB_nz_id = 0; TLB_nz_id < nnz_of_TLB; TLB_nz_id++)
        {
            // 获取特定的bit
            // 当前线程处理第一个非零元的sum_begin_bit，因为右移是带符号的，所以要截一下
            bool cur_sum_begin_bit = (combine_meta >> (31 - TLB_nz_id)) & 0x1;

            // 如果当前非零元是第一个行首非零元，那么之前的结果要存起来，用来跨线程归约的行其他部分中间结果
            if (cur_sum_begin_bit == true)
            {
                if (through_sum_begin_bit == true)
                {
                    // 直接写全局内存，并且如果WLB内TLB的id是0，并且写的位置和TLB的第一个行偏移一致，那就要原子加
                    // y_offset与要写的位置相关
                    if (tid_in_warp == 0 && y_offset == TLB_first_reduce_row)
                    {
                        // 原子加
                        atomicAdd(&(y[WLB_first_row + y_offset]), sum_tmp);
                    }
                    else
                    {
                        y[WLB_first_row + y_offset] = sum_tmp;
                    }
                }
                else
                {
                    // 第一个行首非零元，行其他部分的中间结果
                    row_other_tmp_result = sum_tmp;
                }
            }

            // TLB中间结果写全局内存的行偏移量，这里的行偏移量是相对索引
            y_offset = y_offset + (through_sum_begin_bit & cur_sum_begin_bit);

            // 记录是不是已经经过了sum begin的点
            through_sum_begin_bit = through_sum_begin_bit || cur_sum_begin_bit;
            
            // 根据实际情况清空归约结果
            sum_tmp = cur_sum_begin_bit ? 0 : sum_tmp;

            sum_tmp = sum_tmp + val[global_nz_index] * x[col_index[global_nz_index]];
            
            // 非零元数量
            global_nz_index = global_nz_index + 32;
        }

        // 之后执行线程间归约，首先初始化所有线程的行身中间结果
        // 对于出现过加和起始位置的，那么行身结果应该已经被赋值了，对于没有出现过加和起始位置
        row_other_tmp_result = through_sum_begin_bit ? row_other_tmp_result : sum_tmp;

        // 当前的sum中间值是脑袋，对于没有加和起始点的元素来说脑袋和身体的值是一样的。对于这种情况来说，这个归约是不会被执行的
        row_head_tmp_result = sum_tmp;

        // 也可以用__shfl达成相同的逻辑，不需要为最后一个非零元赋值
        row_other_tmp_result = __shfl_down_sync(0xFFFFFFFF, row_other_tmp_result, 1);
        
        // 最后一个赋值为0，其他的保持不变
        row_other_tmp_result = (tid_in_warp == 31) ? 0 : row_other_tmp_result;

        // 执行一个扫描加，warp最后一个线程的归约偏移肯定是0
        assert(reduce_offset <= 31 - tid_in_warp);

        // 扫描加，scan_tmp_sum中放着扫描加的结果
        float scan_tmp_sum = row_other_tmp_result;
        // 先计算每个线程更低的所有线程row_other_tmp_result的
        float scan_tmp_sum_from_other_thread = __shfl_up_sync(0xFFFFFFFF, scan_tmp_sum, 1);
        scan_tmp_sum = (tid_in_warp >= 1) ? scan_tmp_sum_from_other_thread + scan_tmp_sum : scan_tmp_sum;
        scan_tmp_sum_from_other_thread = __shfl_up_sync(0xFFFFFFFF, scan_tmp_sum, 2);
        scan_tmp_sum = (tid_in_warp >= 2) ? scan_tmp_sum_from_other_thread + scan_tmp_sum : scan_tmp_sum;
        scan_tmp_sum_from_other_thread = __shfl_up_sync(0xFFFFFFFF, scan_tmp_sum, 4);
        scan_tmp_sum = (tid_in_warp >= 4) ? scan_tmp_sum_from_other_thread + scan_tmp_sum : scan_tmp_sum;
        scan_tmp_sum_from_other_thread = __shfl_up_sync(0xFFFFFFFF, scan_tmp_sum, 8);
        scan_tmp_sum = (tid_in_warp >= 8) ? scan_tmp_sum_from_other_thread + scan_tmp_sum : scan_tmp_sum;
        scan_tmp_sum_from_other_thread = __shfl_up_sync(0xFFFFFFFF, scan_tmp_sum, 16);
        scan_tmp_sum = (tid_in_warp >= 16) ? scan_tmp_sum_from_other_thread + scan_tmp_sum : scan_tmp_sum;

        // scan_tmp_sum是当前原子加之后的结果，根据归约偏移量计算归约结果
        scan_tmp_sum_from_other_thread = __shfl_down_sync(0xFFFFFFFF, scan_tmp_sum, reduce_offset);
        row_other_tmp_result = scan_tmp_sum_from_other_thread - scan_tmp_sum + row_other_tmp_result;

        // 对所有线程执行行脑袋和行身体的归约，但是实际上只有带加和起始位置的线程的计算是有效的
        row_head_tmp_result = through_sum_begin_bit ? row_head_tmp_result + row_other_tmp_result : 0;

        // 对于有起始加标记的行来说，才执行写会
        if (through_sum_begin_bit == true)
        {
            // 直接写全局内存，并且如果WLB内TLB的id是0，并且写的位置和TLB的第一个行偏移一致，那就要原子加
            // y_offset与要写的位置相关。如果写的是y的最后一个位置，就要用原子加。最后一个位置需要TLB号和归约偏移是不是加起来是31。
            if ((tid_in_warp == 0 && y_offset == TLB_first_reduce_row) || (tid_in_warp + reduce_offset == 31))
            {
                // 原子加
                atomicAdd(&(y[WLB_first_row + y_offset]), row_head_tmp_result);
            }
            else
            {
                y[WLB_first_row + y_offset] = row_head_tmp_result;
            }
        }
    }
}

// kernal
// 包含两个原始的基本数据
__global__ void kernal(unsigned int* col_index, float* val, float* x, float* y, unsigned int* WLB_row_offset, unsigned int* combine_meta_vec, int bit_flag_bit_num, int reduce_offset_bit_num, int TLB_relative_row_index_bit_num, int WLB_num, int TLB_num, int nnz_of_TLB)
{
    // 查看当前warp号，然后是warp内的线程号
    int global_tid = blockDim.x * blockIdx.x + threadIdx.x;
    int warp_id = global_tid / 32;
    // thread的线程内号
    int tid_in_warp = global_tid % 32;
    int warp_num = blockDim.x * gridDim.x / 32;

    __shared__ unsigned int s_WLB_row_offset[17];
    
    int BLB_id = blockIdx.x;
    int BLB_first_WLB_id = BLB_id * 16;
    
    if (threadIdx.x < 17 && BLB_first_WLB_id + threadIdx.x < WLB_num)
    {
        s_WLB_row_offset[threadIdx.x] = WLB_row_offset[BLB_first_WLB_id + threadIdx.x];
    }

    __syncthreads();
    

    // 处理所有的WLB
    unsigned int WLB_id = warp_id;
    if (WLB_id < WLB_num)
    // for (unsigned int WLB_id = warp_id; WLB_id < WLB_num; WLB_id = WLB_id + warp_num)
    {
        // 当前WLB的首行索引
        unsigned int WLB_first_row = s_WLB_row_offset[WLB_id % 16];
        unsigned int next_WLB_first_row = s_WLB_row_offset[WLB_id % 16 + 1];

        // 对于一个warp一整行的情况，不需要那么复杂，每个TLB执行自己的，并且
        if (WLB_first_row == next_WLB_first_row)
        {
            // 执行简单的整行归约
            unsigned int global_nz_index = WLB_id * nnz_of_TLB * 32 + tid_in_warp;
            // 每个TLB执行加和，然后树状归约
            float tmp_sum = 0;

            for (unsigned int TLB_nz_id = 0; TLB_nz_id < nnz_of_TLB; TLB_nz_id++)
            {
                tmp_sum = tmp_sum + val[global_nz_index] * x[col_index[global_nz_index]];
                global_nz_index = global_nz_index + 32;
            }

            // 执行一个warp内的归约
            #pragma unroll
            for (int offset = 32 >> 1; offset > 0; offset >>= 1)
            {
                tmp_sum += __shfl_xor_sync(0xFFFFFFFF, tmp_sum, offset);
            }
            
            // 第一个线程向共享内存中写内容
            if (tid_in_warp == 0)
            {
                atomicAdd(&(y[WLB_first_row]), tmp_sum);
            }
        }
        else
        {
            // 执行复杂的归约过程
            // 获取合并之后元数据，从这里计算出总TLB号，每个WLB有32个TLB。
            unsigned int global_TLB_id = WLB_id * 32 + tid_in_warp;

            // 获得每一个TLB的元数据
            unsigned int combine_meta = combine_meta_vec[global_TLB_id];

            // 获取每个TLB第一个加和起始位置的相对行索引以及归约偏移，这两个数据都是先右移再左移的方式带来，也可以直接用与运算来处理。
            unsigned int TLB_first_reduce_row = combine_meta << (32 - TLB_relative_row_index_bit_num) >> ((32 - TLB_relative_row_index_bit_num));
            
            // TLB向全局内存写中间结果的行偏移量
            unsigned int y_offset = TLB_first_reduce_row;

            // 获取当前线程的归约偏移量
            unsigned int reduce_offset = combine_meta << (32 - TLB_relative_row_index_bit_num - reduce_offset_bit_num) >> (32 - reduce_offset_bit_num);

            // 用一个变量记录之前是否出现过sum起始点
            bool through_sum_begin_bit = false;

            // 用一个变量存储行脑袋和行身体的中间结果
            float row_head_tmp_result = 0;
            // 除了头之外其他部分的中间结果
            float row_other_tmp_result = 0;
            // 加和的中间结果
            float sum_tmp = 0;

            // 当前非零元的全局索引
            // assert(global_TLB_id < 160);
            unsigned int global_nz_index = WLB_id * nnz_of_TLB * 32 + tid_in_warp;
            
            // 遍历线程所有的非零元
            for (unsigned int TLB_nz_id = 0; TLB_nz_id < nnz_of_TLB; TLB_nz_id++)
            {
                // 获取特定的bit
                // 当前线程处理第一个非零元的sum_begin_bit，因为右移是带符号的，所以要截一下
                bool cur_sum_begin_bit = (combine_meta >> (31 - TLB_nz_id)) & 0x1;

                // 如果当前非零元是第一个行首非零元，那么之前的结果要存起来，用来跨线程归约的行其他部分中间结果
                if (cur_sum_begin_bit == true)
                {
                    if (through_sum_begin_bit == true)
                    {
                        // 直接写全局内存，并且如果WLB内TLB的id是0，并且写的位置和TLB的第一个行偏移一致，那就要原子加
                        // y_offset与要写的位置相关
                        if (tid_in_warp == 0 && y_offset == TLB_first_reduce_row)
                        {
                            // 原子加
                            atomicAdd(&(y[WLB_first_row + y_offset]), sum_tmp);
                        }
                        else
                        {
                            y[WLB_first_row + y_offset] = sum_tmp;
                        }
                    }
                    else
                    {
                        // 第一个行首非零元，行其他部分的中间结果
                        row_other_tmp_result = sum_tmp;
                    }
                }

                // TLB中间结果写全局内存的行偏移量，这里的行偏移量是相对索引
                y_offset = y_offset + (through_sum_begin_bit & cur_sum_begin_bit);

                // 记录是不是已经经过了sum begin的点
                through_sum_begin_bit = through_sum_begin_bit || cur_sum_begin_bit;
                
                // 根据实际情况清空归约结果
                sum_tmp = cur_sum_begin_bit ? 0 : sum_tmp;

                sum_tmp = sum_tmp + val[global_nz_index] * x[col_index[global_nz_index]];
                
                // 非零元数量
                global_nz_index = global_nz_index + 32;
            }

            // 之后执行线程间归约，首先初始化所有线程的行身中间结果
            // 对于出现过加和起始位置的，那么行身结果应该已经被赋值了，对于没有出现过加和起始位置
            row_other_tmp_result = through_sum_begin_bit ? row_other_tmp_result : sum_tmp;

            // 当前的sum中间值是脑袋，对于没有加和起始点的元素来说脑袋和身体的值是一样的。对于这种情况来说，这个归约是不会被执行的
            row_head_tmp_result = sum_tmp;

            // 也可以用__shfl达成相同的逻辑，不需要为最后一个非零元赋值
            row_other_tmp_result = __shfl_down_sync(0xFFFFFFFF, row_other_tmp_result, 1);
            
            // 最后一个赋值为0，其他的保持不变
            row_other_tmp_result = (tid_in_warp == 31) ? 0 : row_other_tmp_result;

            // 执行一个扫描加，warp最后一个线程的归约偏移肯定是0
            assert(reduce_offset <= 31 - tid_in_warp);

            // 扫描加，scan_tmp_sum中放着扫描加的结果
            float scan_tmp_sum = row_other_tmp_result;
            // 先计算每个线程更低的所有线程row_other_tmp_result的
            float scan_tmp_sum_from_other_thread = __shfl_up_sync(0xFFFFFFFF, scan_tmp_sum, 1);
            scan_tmp_sum = (tid_in_warp >= 1) ? scan_tmp_sum_from_other_thread + scan_tmp_sum : scan_tmp_sum;
            scan_tmp_sum_from_other_thread = __shfl_up_sync(0xFFFFFFFF, scan_tmp_sum, 2);
            scan_tmp_sum = (tid_in_warp >= 2) ? scan_tmp_sum_from_other_thread + scan_tmp_sum : scan_tmp_sum;
            scan_tmp_sum_from_other_thread = __shfl_up_sync(0xFFFFFFFF, scan_tmp_sum, 4);
            scan_tmp_sum = (tid_in_warp >= 4) ? scan_tmp_sum_from_other_thread + scan_tmp_sum : scan_tmp_sum;
            scan_tmp_sum_from_other_thread = __shfl_up_sync(0xFFFFFFFF, scan_tmp_sum, 8);
            scan_tmp_sum = (tid_in_warp >= 8) ? scan_tmp_sum_from_other_thread + scan_tmp_sum : scan_tmp_sum;
            scan_tmp_sum_from_other_thread = __shfl_up_sync(0xFFFFFFFF, scan_tmp_sum, 16);
            scan_tmp_sum = (tid_in_warp >= 16) ? scan_tmp_sum_from_other_thread + scan_tmp_sum : scan_tmp_sum;

            // scan_tmp_sum是当前原子加之后的结果，根据归约偏移量计算归约结果
            scan_tmp_sum_from_other_thread = __shfl_down_sync(0xFFFFFFFF, scan_tmp_sum, reduce_offset);
            row_other_tmp_result = scan_tmp_sum_from_other_thread - scan_tmp_sum + row_other_tmp_result;

            // 对所有线程执行行脑袋和行身体的归约，但是实际上只有带加和起始位置的线程的计算是有效的
            row_head_tmp_result = through_sum_begin_bit ? row_head_tmp_result + row_other_tmp_result : 0;

            // 对于有起始加标记的行来说，才执行写会
            if (through_sum_begin_bit == true)
            {
                // 直接写全局内存，并且如果WLB内TLB的id是0，并且写的位置和TLB的第一个行偏移一致，那就要原子加
                // y_offset与要写的位置相关。如果写的是y的最后一个位置，就要用原子加。最后一个位置需要TLB号和归约偏移是不是加起来是31。
                if ((tid_in_warp == 0 && y_offset == TLB_first_reduce_row) || (tid_in_warp + reduce_offset == 31))
                {
                    // 原子加
                    atomicAdd(&(y[WLB_first_row + y_offset]), row_head_tmp_result);
                }
                else
                {
                    y[WLB_first_row + y_offset] = row_head_tmp_result;
                }
            }
        }
    }
}

int main()
{
    // 做一个一共n行的矩阵，每一行的非零元数量为n，列号一个个增加，值数组全为1
    // 每个线程的非零元数量
    int nnz_of_each_thread = 5;

    vector<unsigned int> row_index;
    vector<unsigned int> col_index;
    vector<float> val_arr;

    unsigned int col_num = 0;
    unsigned int row_num = 0;
    unsigned int nnz = 0;

    // get_coo_vector_of_matrix_from_file(string file_name, unsigned int* col_num, unsigned int* row_num, unsigned int* nnz, vector<unsigned int> &row_index_arr, vector<unsigned int> &col_index_arr, vector<float> &val_arr)
    get_coo_vector_of_matrix_from_file("/home/duzhen/spmv_csr5_2021/resource/webbase-1M.mtx.coo", &col_num, &row_num, &nnz, row_index, col_index, val_arr);

    cout << nnz << endl;
    cout << row_num << endl;
    cout << col_num << endl;

    if (nnz % (32 * nnz_of_each_thread) != 0)
    {
        // 换成一个新的值
        nnz = (nnz / (32 * nnz_of_each_thread) + 1) * 32 * nnz_of_each_thread;
        // 执行一个padding
        // padding_coo_matrix_to_target_nnz(vector<unsigned int> &row_index_arr, vector<unsigned int> &col_index_arr, vector<float> &val_arr, unsigned int target_nnz)
        padding_coo_matrix_to_target_nnz(row_index, col_index, val_arr, nnz);
    }

    cout << nnz << endl;

    // x向量
    vector<float> x;
    vector<float> y;

        
    
    // // 总非零元数量5*32
    // int warp_num = 5;
    // int nnz = 32 * warp_num * nnz_of_each_thread;
    
    // int cur_nnz = 0;
    // for (int row_id = 0; cur_nnz < nnz; row_id++)
    // {
    //     // 给列号赋值
    //     // for (int col_id = 0; col_id < row_id + 1; col_id++)
    //     // {
    //     //     row_index.push_back(row_id);
    //     //     col_index.push_back(col_id);
    //     //     val_arr.push_back(1);
    //     //     cur_nnz++;
            
    //     //     if (cur_nnz == nnz)
    //     //     {
    //     //         break;
    //     //     }
    //     // }

    //     // 每一行列的数量相同
    //     for (int col_id = 0; col_id < 1; col_id++)
    //     {
    //         row_index.push_back(row_id);
    //         col_index.push_back(col_id);
    //         val_arr.push_back(1);

    //         cur_nnz++;

    //         if (cur_nnz == nnz)
    //         {
    //             break;
    //         }
    //     }
    // }

    // 初始化结果
    for (int row_id = 0; row_id <= row_num; row_id++)
    {
        y.push_back(0);
    }

    // 初始化向量
    for (int col_id = 0; col_id <= col_num; col_id++)
    {
        x.push_back(1);
    }

    // // 再添加一些
    // for (int i = 0; i < 100; i++)
    // {
    //     x.push_back(1);
    // }

    // assert(row_index.size() % 32 == 0);

    // // cout << row_index.size() << endl;
    // // 打印COO格式的矩阵
    // // print_coo_matrix(row_index, col_index, val_arr);
    // store_coo_matrix_to_file(row_index, col_index, val_arr, "/home/duzhen/spmv_csr5_2021/coo_matrix");
    
    // // 打印标记
    vector<vector<bool>> bool_flag_of_sum_begin = get_sum_begin_bool_flag_of_each_thread(row_index, nnz_of_each_thread);

    // store_bool_flag_of_sum_begin_to_file(bool_flag_of_sum_begin, "/home/duzhen/spmv_csr5_2021/bool_flag_of_sum_begin");

    // // WLB首行
    vector<unsigned int> WLB_row_offset = first_global_row_index_of_each_warp(row_index, nnz_of_each_thread);

    // store_vec_to_file<unsigned int>(WLB_row_offset, "/home/duzhen/spmv_csr5_2021/WLB_row_offset");

    vector<unsigned int> TLB_row_offset = first_relative_row_index_of_each_thread(row_index, WLB_row_offset, bool_flag_of_sum_begin, nnz_of_each_thread);

    // store_vec_to_file<unsigned int>(TLB_row_offset, "/home/duzhen/spmv_csr5_2021/TLB_row_offset");

    vector<unsigned int> reduce_offset_vec = get_reduce_offset_vec(bool_flag_of_sum_begin);
    
    // store_vec_to_file<unsigned int>(reduce_offset_vec, "/home/duzhen/spmv_csr5_2021/reduce_offset_vec");

    int bit_num_of_bit_flag;
    int bit_num_of_reduce_offset;
    int bit_num_of_relative_row_index;

    vector<unsigned int> combine_meta_vec = combine_meta_to_a_unsigned_int_vec(bool_flag_of_sum_begin, reduce_offset_vec, TLB_row_offset, &bit_num_of_bit_flag, &bit_num_of_reduce_offset, &bit_num_of_relative_row_index);

    // store_combine_meta_vec(combine_meta_vec, "/home/duzhen/spmv_csr5_2021/combine_meta_vec");

    col_index = stagger_col_index_in_warp(col_index);
    val_arr = stagger_val_in_warp(val_arr);

    // // 几个数据
    // // unsigned int* col_index, float* val, unsigned int* WLB_row_offset, unsigned int* combine_meta_vec
    unsigned int *dcol_index = NULL;
    float *dval_arr = NULL;
    unsigned int *dcombine_meta_vec = NULL;
    unsigned int *dWLB_row_offset = NULL;
    float *dx = NULL;
    float* dy = NULL;

    cudaMalloc(&dcol_index, sizeof(unsigned int) * col_index.size());
    cudaMalloc(&dval_arr, sizeof(float) * val_arr.size());
    cudaMalloc(&dcombine_meta_vec, sizeof(unsigned int) * combine_meta_vec.size());
    cudaMalloc(&dWLB_row_offset, sizeof(unsigned int) * WLB_row_offset.size());
    cudaMalloc(&dx, sizeof(float) * x.size());
    cudaMalloc(&dy, sizeof(float) * y.size());

    // // 执行一个拷贝
    cudaMemcpy(dcol_index, &(col_index[0]), sizeof(unsigned int) * col_index.size(), cudaMemcpyHostToDevice);
    cudaMemcpy(dval_arr, &(val_arr[0]), sizeof(float) * val_arr.size(), cudaMemcpyHostToDevice);
    cudaMemcpy(dcombine_meta_vec, &(combine_meta_vec[0]), sizeof(unsigned int) * combine_meta_vec.size(), cudaMemcpyHostToDevice);
    cudaMemcpy(dWLB_row_offset, &(WLB_row_offset[0]), sizeof(unsigned int) * WLB_row_offset.size(), cudaMemcpyHostToDevice);
    cudaMemcpy(dx, &(x[0]), sizeof(float) * x.size(), cudaMemcpyHostToDevice);
    cudaMemcpy(dy, &(y[0]), sizeof(float) * y.size(), cudaMemcpyHostToDevice);

    cout << "TLB_num:" << combine_meta_vec.size() << endl;
    cout << "nnz:" << col_index.size() << endl;
    cudaDeviceSynchronize();

    unsigned int TLB_num = combine_meta_vec.size();
    assert(TLB_num % 32 == 0);
    
    // 实际活跃的线程会稍微多一点点
    unsigned int BLB_num;

    if (TLB_num % 512 == 0)
    {
        BLB_num = TLB_num / 512;
    }
    else
    {
        BLB_num = TLB_num / 512 + 1;
    }
    
    unsigned int repeat_num = 32000;

    struct timeval start,end;
    gettimeofday(&start, NULL);

    // 遍历特定的数量，计算gflops
    for (int repeat_id = 0; repeat_id < repeat_num; repeat_id++)
    {
        kernal_basic<<<BLB_num, 512>>>(dcol_index, dval_arr, dx, dy, dWLB_row_offset, dcombine_meta_vec, bit_num_of_bit_flag, bit_num_of_reduce_offset, bit_num_of_relative_row_index, WLB_row_offset.size(), combine_meta_vec.size(), nnz_of_each_thread);
        cudaDeviceSynchronize();
    }
    // // __global__ void kernal(unsigned int* col_index, float* val, float* x, float* y, unsigned int* WLB_row_offset, unsigned int* combine_meta_vec, int bit_flag_bit_num, int reduce_offset_bit_num, int TLB_relative_row_index_bit_num, int WLB_num, int TLB_num, int nnz_of_TLB)
    gettimeofday(&end, NULL);

    long timeuse = 1000000 * (end.tv_sec - start.tv_sec ) + end.tv_usec - start.tv_usec;
    double gflops = ((double)2.0 * nnz * repeat_num / ((double)timeuse / 1000000)) / 1000000000;

    printf("time=%fms, gflops=%f\n",timeuse /1000.0, gflops);

    cudaError_t cuda_err = cudaGetLastError();

    if(cudaSuccess != cuda_err)
    {
        fprintf(stderr,"%s\n", cudaGetErrorString(cuda_err));

        return -1;
    }

    // // 将y拷贝出来
    cudaMemcpy(&(y[0]), dy, sizeof(float) * y.size(), cudaMemcpyDeviceToHost);

    // // 将y的内容输出到文件中
    store_vec_to_file<float>(y, "/home/duzhen/spmv_csr5_2021/tmp/y");

    // // 
}
